Meshfree Algorithm for Level Set Evolution

ABSTRACT

The present invention is a system and method for simulating the motion of an interface. The interface moving through a simulation space. The invention includes simulating the interface using a level set function to describe a position and shape of the interface in the simulation space at a first point in time. The invention also includes describing the level set function at the first point time using a meshfree method. The invention further includes describing a motion of the interface from the first point in time to a second point time using a level set evolution method. The invention also includes finding an approximate solution to the level set evolution method using the meshfree method to describe the level set function at the second point in time.

CROSS-REFERENCE TO RELATED APPLICATIONS

U.S. patent application Ser. No. 11/031,906 (Attorney Docket No.AP210HO) filed on Jan. 6, 2005, is hereby incorporated by reference inits entirety.

BACKGROUND

1. Field of the Invention

The present invention relates generally to systems and methods forevaluating level set functions and their evolution over time. Whereinthe level set function represents a physical phenomenon.

2. Description of the Related Art

The level set method is a robust method for accurately describing anarbitrarily complex interface. The level set approach defines a smoothsigned distance function φ({right arrow over (x)},t), which representsthe interface as a zero level set, S=φ({right arrow over (x)},t)=0.Instead of explicitly tracking the interface, the interface is now being“captured” by calculating the interface for which φ({right arrow over(x)},t) is equal to zero. Generalizing to higher dimensions. The zerolevel set S is an n dimensional manifold in an n+1 dimensional space.The level set function φ is defined throughout the n+1 dimensionalspace. The absolute value of the level set function φ at any point inspace is equal to the shortest distance between that point in space andthe zero level set S. The sign of the level set function φ is positiveon one side of the zero level set S and negative on the other side ofthe zero level set S.

An advantage of using the level set method over other methods is thelevel set method's ability to handle changes in the topology such asinterface break-ups and merges. The level set method has demonstratedits ability to help the study of complex multi-fluid systems and variousfluid dynamics problems such as the motion of droplets and bubbles, freesurface flow, inflation of an airbag, etc. and multidimensional imageprocessing problems such as surface processing and surfacereconstruction.

The finite element (FE) method is a common method for solving a partialdifferential equation such as a level set evolution equation. The levelset evolution equation describes changes that occur in the level setfunction φ over an additional dimension such as time. Finite differencemethods and spectral methods may also be used to solve the level setevolution equation. The FE method solves the level set evolutionequation by capturing the interface S and by interpolating finiteelement shape functions, much as other state variables such as pressure,velocity, or temperature are interpolated.

The interface capturing technique based on the FE method can be robustand simple to implement. However, it requires high resolution tooeffectively capture the evolution of the interface. The FE methodinvolves the discretization of a given problem domain into simplegeometric shapes called elements. The accuracy of the FE method dependson the size of the elements. Calculating the evolution of a level setfunction using the FE method with sufficient accuracy requires solvinglarge systems of linear equations. Solving large systems of linearequations can be computationally expensive.

The present invention is aimed at reducing the computational resourcesrequired to solve level set evolution equations relative to the finiteelement method.

SUMMARY OF THE INVENTION

The present invention is a system and method for simulating the motionof an interface. The interface is moving through a simulation space. Theinvention may include simulating the interface using a level setfunction to describe a position and shape of the interface in thesimulation space at a first point in time. The invention may alsoinclude describing the level set function at the first point time usinga meshfree method. The invention may further include describing a motionof the interface from the first point in time to a second point time,using a level set evolution method. The invention may also includefinding an approximate solution to the level set evolution method usingthe meshfree method to describe the level set function at the secondpoint in time.

In an embodiment of the present invention, the meshfree method may be aReproducing Kernel Particle Method. In an embodiment of the presentinvention, the interface may be between a first fluid and a secondfluid. In an embodiment of the present invention, the first fluid may beink, and the second fluid may be air.

An embodiment of the present invention may be a computer-readable mediumencoded with instructions for simulating the motion of an interface. Anembodiment of the present invention may include a system including amemory module and instructions for simulating the motion of aninterface.

In an embodiment of the present invention, the meshfree method mayapproximate the level set function using a first family of functionsthat define the smoothness of the approximation of the level setfunction. Each function in the first family of functions may have acompact support. The size of the compact support may be a function of amaximum of a set of k-th ordered statistics of a plurality of distancesbetween sets of distinct nodes in the simulation space. In an embodimentof the present invention, k may be selected from the group consisting of2, 3, 4, and 5. In an embodiment of the present invention, the size ofthe compact support may also be a function of an order of an enrichmentfunction. In an embodiment of the present invention, each function inthe family of function may be positive in the area bounded by the sizeof the compact support.

In an embodiment of the present invention, the level set evolutionmethod may be described in terms of spatial derivates of the level setfunction. The spatial derivatives of the level set function may bediscretized. A mesh free method may be used to describe each of thediscretized terms of the spatial derivatives of the level set function.

In an embodiment of the present invention, the level set evolutionmethod may include an integral along a boundary surrounding thesimulation space. The integral may include a spatial derivate of thelevel set function, a test function, and components of a vector normalto the boundary.

An embodiment of the present invention may include designing a physicalobject based on the results of the method of simulating the motion of aninterface. In an embodiment of the present invention, the physicalobject may be an ink jet head. An embodiment of the present inventionmay be the physical object that is designed based on the results of themethod of simulating the motion of an interface.

An embodiment of the present may include a meshfree method thatapproximates the level set function using a first family of functionsthat define the smoothness of the approximation of the level setfunction. Each function in the first family of functions may have acompact support The size of the compact support may be a function of amaximum of a set of k-th ordered statistics of a plurality of power sumsof differences between components that define the positions of sets ofdistinct nodes in the simulation space.

Other objects and attainments together with a fuller understanding ofthe invention will become apparent and appreciated by referring to thefollowing description and claims taken in conjunction with theaccompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

In the drawings wherein like reference symbols refer to like parts.

FIG. 1A is an illustration of a FE approximation of a sine function anda derivative of the FE approximation;

FIG. 1B is an illustration of a Meshfree approximation of a sinefunction and a derivate of the Meshfree approximation of the sinefunction;

FIGS. 2A-2C show the evolution of a one-dimensional level set functionover a series of time steps;

FIGS. 3A-3E are illustrations of the results of an embodiment of thepresent invention used to simulate the evolution of a two-dimensionallevel set function;

FIGS. 4A-B are illustrations of the evolution of one dimensional levelset function in which the level set is not reinitialized and theboundary integral is ignored;

FIGS. 5A-B are illustrations of the evolution of one dimensional levelset function in which the level set is reinitialized and the boundaryintegral is ignored;

FIGS. 6A-F are illustration of the evolution of a two dimensional levelset function φ at several different time steps in which the level setfunction is reinitialized as used in an embodiment of the presentinvention;

FIG. 7 is an illustration of a boundary;

FIGS. 8A-B are illustrations of the evolution of one dimensional levelset function in which the level set is reinitialized and the boundaryintegral is not ignored; and

FIG. 9 is an illustration of a system in which an embodiment of thesystem may be practiced.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

The present invention minimizes the computing cost required to solve alevel set evolution equation. In the present invention, a Meshfreemethod is used to discretize a problem domain. Instead of discretizingthe problem domain as in the FE method, the meshfree method scattersnodes throughout the problem domain. The Meshfree method allows a globallevel of approximation and eliminates the use of elements.

The Meshfree method has several advantages over the FE method. Among theadvantages of using, the Meshfree method over the FE method is thatsmooth derivatives of the shape function may be obtained over the entirenumerical domain. FIG. 1A is an illustration of a FE approximation of asine function and a derivative of the FE approximation. FIG. 1B is anillustration of a Meshfree approximation of a sine function and aderivate of the Meshfree approximation. As seen in FIG. 1A the derivateof the FE method approximation includes multiple discontinuities.Whereas the derivative of the Meshfree approximation as shown in FIG. 1Bis smooth.

Galerkin methods typically include the calculation of both thederivative of the shape function and the shape function itself. Thus,Galerkin methods benefit when the derivate is well behaved.

Meshfree Method

An embodiment of the present invention makes use of a Reproducing KernelParticle Method (RPKM) much like the method disclosed in U.S. patentapplication Ser. No. 11/031,906 (AP210HO) filed on Jan. 6, 2005 which ishereby incorporated by reference in it's entirety. An embodiment of thepresent invention may make use of RKPM with monomial basis functions.

A variable u({right arrow over (x)}) may be approximated with a discretereproducing kernel approximation, denoted by u^(h), as shown in equation(1).

$\begin{matrix}{{u^{h}\left( \overset{->}{x} \right)} = {\sum\limits_{I = 1}^{NP}{{{\overset{\_}{\Phi}}_{a}\left( {\overset{->}{x};{\overset{->}{x} - {\overset{->}{x}}_{I}}} \right)}d_{I}}}} & (1)\end{matrix}$

In equation (1), NP is the number of discrete nodes that have beenscattered throughout the simulation space. The variable d₁ representsthe coefficients of the approximation. The function Φ _(a)({right arrowover (x)};{right arrow over (x)}−{right arrow over (x)}_(i)) is thereproducing kernel and is formed by a multiplication of two functions asshown in equation (2).

Φ({right arrow over (x)};{right arrow over (x)}−{right arrow over (x)}_(I))=C({right arrow over (x)};{right arrow over (x)}−{right arrow over(x)} _(I))Φ_(a)({right arrow over (x)}−{right arrow over (x)} _(I))  (2)

As used in equation (2) the function Φ_(a)({right arrow over (x)}−{rightarrow over (x)}_(i)) is a kernel function that defines the smoothness ofthe approximation and has a compact support as measured by thecoefficient a. The function C({right arrow over (x)};{right arrow over(x)}−{right arrow over (x)}_(i)) acts as an enrichment function (i.e., acorrection function). The enrichment function is used to satisfy then-th order reproducing conditions as stated in equation (3). Thevariable x_(iI) is the nodal value of x_(i) at node I.

$\begin{matrix}{{{\sum\limits_{I = 1}^{NP}{{\overset{\_}{\Phi}\left( {\overset{->}{x};{\overset{->}{x} - {\overset{->}{x}}_{I}}} \right)}x_{1\; I}^{p}x_{2\; I}^{q}x_{3\; I}^{r}}} = {x_{1}^{p}x_{2}^{q}x_{3}^{q}}}{{{p + q + r} = {0\mspace{14mu} \ldots \mspace{14mu} n}};\left( {{x_{1} \equiv x},{x_{2} \equiv y},{x_{3} \equiv z}} \right)}} & (3)\end{matrix}$

The n-th order reproducing function described in equation 3 has beendefined in terms of the three dimensions: x₁; x₂; and x₃, an individualskilled in the art would appreciate how to modify equation (3) for bothhigher and lower dimensions. To meet the n-th order reproducingconditions of equation (3), the enrichment function C({right arrow over(x)};{right arrow over (x)}−{right arrow over (x)}_(I)) may beconstructed as a linear combination of complete n-th order monomialfunctions as shown in equation (4).

$\begin{matrix}\begin{matrix}{{C\left( {\overset{->}{x};{\overset{->}{x} - {\overset{->}{x}}_{I}}} \right)} = {\sum\limits_{{p + q + r} = 0}^{n}{\left( {x_{1} - x_{1I}} \right)^{p}\left( {x_{2} - x_{2\; I}} \right)^{q}\left( {x_{3} - x_{3\; I}} \right)^{r}{b_{pqr}\left( \overset{->}{x} \right)}}}} \\{\equiv {{H^{T}\left( {\overset{->}{x} - {\overset{->}{x}}_{I}} \right)}{b\left( \overset{->}{x} \right)}}}\end{matrix} & (4)\end{matrix}$

The functions b_(pqr)({right arrow over (x)}) are the coefficients ofthe monomial basis functions and are functions of the spatial vector{right arrow over (x)}, The matrix b({right arrow over (x)}) is an m×1matrix of the functions b_(pqr)({right arrow over (x)}). Wherein m isthe number of elements in the polynomial listed in equation (4). Thematrix H({right arrow over (x)}−{right arrow over (x)}_(I)) is an m×1matrix containing the monomial basis functions as shown in equation (5).

H ^(T)({right arrow over (x)}−{right arrow over (x)} _(I))=[1, x ₁ −x_(1I) , x ₂ −x _(2I) , x ₃ −x _(3I), (x ₁ −x _(1I))², . . . , (x ₃ −x_(3I))^(n)]  (5)

Using the notation above, n-th order reproducing conditions, equation(3) can be rewritten as equation (6).

$\begin{matrix}{{\sum\limits_{I = 1}^{NP}{{H\left( {\overset{->}{x} - {\overset{->}{x}}_{I}} \right)}{\overset{\_}{\Phi}\left( {\overset{->}{x};{\overset{->}{x} - {\overset{->}{x}}_{I}}} \right)}}} = {{H(0)} = \left\lbrack {1,0,\ldots \mspace{14mu},0} \right\rbrack}} & (6)\end{matrix}$

The reproducing conditions of equation (6) may be written as matrixequation (7) which may be used to solve for the coefficients b({rightarrow over (x)}).

M({right arrow over (x)})b({right arrow over (x)})=H(0)   (7)

The moment matrix M({right arrow over (x)}) may be defined in the mannershown in equation (8).

$\begin{matrix}{{M\left( \overset{->}{x} \right)} = {\sum\limits_{I = 1}^{NP}{{H\left( {\overset{->}{x} - {\overset{->}{x}}_{I}} \right)}{H^{T}\left( {\overset{->}{x} - {\overset{->}{x}}_{I}} \right)}{\Phi_{a}\left( {\overset{->}{x} - {\overset{->}{x}}_{I}} \right)}}}} & (8)\end{matrix}$

M({right arrow over (x)}) is the moment matrix of the functionΦ_(a)({right arrow over (x)}−{right arrow over (x)}_(I)). In anembodiment of the present invention the matrix M({right arrow over (x)})may be required to be invertible. In order to help insure that thematrix M({right arrow over (x)}) in equation (8) is invertible, thesupport a of the function φ_(a)({right arrow over (x)}−{right arrow over(x)}_(I)) needs to be greater than a minimum size that is related to theorder m of basis functions used in the enrichment function C({rightarrow over (x)};{right arrow over (x)}−{right arrow over (x)}_(I)) andthe nodal spacing, and Φ_(a)({right arrow over (x)}−{right arrow over(x)}_(I)) must be a positive function within the support.

For example, the support a may be a function of the distance between anytwo distinct nodes; h_(IJ)=|{right arrow over (x)}_(I) −{right arrowover (x)} _(J)|. The matrix D may be defined as complete set of thesedistances. A set of minimum distances between any two nodes may bedefined as minimum of the matrix D across one index of the matrix

${\min\limits_{I}(D)} = {{\min\limits_{I}\left( {{{\overset{\rightarrow}{x}}_{I} - {\overset{\rightarrow}{x}}_{J}}} \right)} = {D_{{(1)}_{I}}.}}$

In an alternative embodiment of the present invention, a k-th orderedstatistic of the matrix D may be used instead of the minimum wherein kis 2, 3, or 4; D_((k)) _(I) . The support a may be a function f of amaximum of the k-th ordered statistic of the matrix

$D;{{\max\limits_{J}\left( D_{{(k)}_{I}} \right)} = \left\lfloor D_{{(k)}_{I}} \right\rfloor_{{({NP})}_{J}}};$

a = f([D_((k)_(I))]_((NP)_(J))), .

An individual skilled in the art will appreciate that all elements ofthe matrix D need not be calculated to obtain the information todetermine the support a An individual skilled in the art will appreciatethat the diagonal elements of the matrix D are zero, and need not beused in the calculation described above. In an alternative embodiment ofthe present invention the function f may also be a function of m theorder of the basis functions, used in the enrichment function

${C\left( {\overset{\rightarrow}{x};{\overset{\rightarrow}{x} - {\overset{\rightarrow}{x}}_{I}}} \right)};{a = {{f\left( {\left\lbrack D_{{(k)}_{I}} \right\rbrack_{{({NP})}_{J}},m} \right)}.}}$

An alternative embodiment of the embodiment of the present invention mayuse a different method to approximate the distance between two nodes.For example it may be shown that this

$a = {f\left( {\sqrt{\left\lbrack \left( {\sum\limits_{i}^{\;}\left( {x_{iI} - x_{iJ}} \right)^{2}} \right)_{{(k)}_{I}} \right\rbrack_{{({NP})}_{J}}},m} \right)}$

method of calculating the support a is equivalent to the methoddescribed above. A further simplification of the calculation may includeremoving the square root;

$a = {{f\left( {\left\lbrack \left( {\sum\limits_{i}^{\;}\left( {x_{iI} - x_{iJ}} \right)^{2}} \right)_{{(k)}_{I}} \right\rbrack_{{({NP})}_{J}},m} \right)}.}$

A generalization of this method may be written in terms of power sums ofdifference

$a = {f\left( {\left\lbrack \left( {\sum\limits_{i}^{\;}\left( {x_{iI} - x_{iJ}} \right)^{2\; n}} \right)_{{(k)}_{I}} \right\rbrack_{{({NP})}_{J}},m} \right)}$

or

$a = {{f\left( {\left\lbrack \left( {\sum\limits_{i}^{\;}{{x_{iI} - x_{iJ}}}^{{2\; n} + 1}} \right)_{{(k)}_{I}} \right\rbrack_{{({NP})}_{J}},m} \right)}.}$

In this generalization, n may be any integer. For example, n may beselected from the group consisting of 1, 2, 3, 4, 5, and 6. As n getslarger greater distances have a larger effect.

Equation (1) can be rewritten as equation (9).

$\begin{matrix}{u^{h} = {\sum\limits_{I = 1}^{NP}{{\Psi_{I}\left( \overset{->}{x} \right)}d_{I}}}} & (9)\end{matrix}$

In which the function Ψ_(I)({right arrow over (x)}) is solved usingequations (2), (4), and (7). Equation (10) is a definition ofΨ_(I)({right arrow over (x)}) in terms of the reproducing kernelapproximation.

Ψ_(I)({right arrow over (x)})=H ^(T)(0)M ⁻¹({right arrow over(x)})H({right arrow over (x)}−{right arrow over (x)} _(I))Φ_(a)({rightarrow over (x)}−{right arrow over (x)}_(I))   (10)

When monomial basis functions are used in the reproducing kernelapproximation, the smoothness and compact support properties of theshape function Ψ_(I)({right arrow over (x)}) are identical to those ofthe kernel function Φ_(a)({right arrow over (x)}−{right arrow over(x)}_(I)). Multi-dimensional kernel functions may be constructed byusing the product of one-dimensional shape functions. Alternatively, themulti-dimensional kernel function may be constructing considering thedistance between nodes |{right arrow over (x)}−{right arrow over(x)}_(I)| as an independent variable in the evaluation of the kernelfunctions.

Level Set Formulation and Temporal Discretization

In the level set method, the free surface is identified as a zero levelset, i.e. φ({right arrow over (x)},t)=0. The free surface moves throughspace and time. A speed function u_(i) can be used to describe how thefree surface moves. Equation (12) is a level set evolution equation thatdescribes how the level set function φ evolves over time in terms of thespeed function u_(i).

$\begin{matrix}{{\frac{\partial\varphi}{\partial t} + {\sum\limits_{i}^{\;}{u_{i}\frac{\partial\varphi}{\partial x_{i}}}}} = 0} & (12)\end{matrix}$

An embodiment of the present invention may be directed toward a methodfor determining the level set function at a future time given thecurrent state of the level set function and equation (12). Prior artmethods have suffered from instabilities when standard discretizationtechniques are used to solve convection dominated problems such as thenone described in equation (12). An embodiment of the present inventionaddresses this issue by making use of the Reproducing Kernel ParticleMethod described above in combination with an explicit Taylor Galerkindiscretization. For example, the level set function φ may be expanded byusing a Taylor series in time. Equation (13) is an example of thisexpansion in which the terms of the second order are retained.

$\begin{matrix}{\varphi^{n + 1} = {\varphi^{n} + {\Delta \; t\frac{\partial\varphi^{n}}{\partial t}} + {\frac{\Delta \; t^{2}}{2}\frac{\partial^{2}\varphi^{n}}{\partial t^{2}}}}} & (13)\end{matrix}$

From equation (12) may be rearranged as equation (14) showing therelationship between the spatial derivative of the level set functionand a temporal derivative of the level set function.

$\begin{matrix}{\frac{\partial\varphi^{n}}{\partial t} = {- {\sum\limits_{i}^{\;}{u_{i}\frac{\partial\varphi^{n}}{\partial x_{i}}}}}} & (14)\end{matrix}$

The second derivate of the level set function is a reasonable extensionof equation (14).

$\begin{matrix}{\frac{\partial^{2}\varphi^{n}}{\partial t^{2}} = {{- {\sum\limits_{i}^{\;}{u_{i}\frac{\partial}{\partial x_{i}}\left( \frac{\partial\varphi^{n}}{\partial t} \right)}}} = {\sum\limits_{i,j}^{\;}{u_{i}\frac{\partial}{\partial x_{i}}\left( {u_{j}\frac{\partial\varphi^{n}}{\partial x_{j}}} \right)}}}} & (15)\end{matrix}$

Using equations (14) and (15) equation (13) can be rewritten solely interms of spatial derivates as shown in equation (16).

$\begin{matrix}{\varphi^{n + 1} = {\varphi^{n} - {\Delta \; t{\sum\limits_{i}^{\;}\left\lbrack {u_{i}\frac{\partial\varphi}{\partial x_{i}}} \right\rbrack^{n}}} + {\frac{\Delta \; t^{2}}{2}{\sum\limits_{i,j}^{\;}\left\lbrack {u_{i}\frac{\partial}{\partial x_{i}}\left( {u_{j}\frac{\partial\varphi}{\partial x_{j}}} \right)} \right\rbrack^{n}}}}} & (16)\end{matrix}$

Equation (16) may than be written in a form that is more amenable toReproducing Kernel Particle Method as shown in equation (17)

$\begin{matrix}\begin{matrix}{{\Delta \; \varphi} = {\varphi^{n + 1} - \varphi^{n}}} \\{= {{{- \Delta}\; t{\sum\limits_{i}^{\;}\left\lbrack {u_{i}\frac{\partial\varphi}{\partial x_{i}}} \right\rbrack^{n}}} + {\frac{\Delta \; t^{2}}{2}{\sum\limits_{i,j}^{\;}\left\lbrack {u_{i}\frac{\partial}{\partial x_{i}}\left( {u_{j}\frac{\partial\varphi}{\partial x_{j}}} \right)} \right\rbrack^{n}}}}}\end{matrix} & (17)\end{matrix}$

The level set function may than be rewritten as in equation (18) interms of a test function Ψ.

φ=Ψ^(T)Φ  (18)

Equation (18) may than be substituted into equation (17) as the whichalso integrated over the solution space in terms of the test function Ψ,as shown in equation (19).

$\begin{matrix}{{\int_{\Omega}^{\;}{\Psi \; \Psi^{T}\ {{\Omega \left( {\Phi^{n + 1} - \Phi^{n}} \right)}}}} = \begin{matrix}{{{- \Delta}\; t{\int_{\Omega}^{\;}{\sum\limits_{i}^{\;}{\Psi \; u_{i}^{n}\frac{\partial\varphi^{n}}{\partial x_{i}}\ {\Omega}}}}} +} \\{\frac{\Delta \; t^{2}}{2}{\int_{\Omega}^{\;}{\sum\limits_{i,j}^{\;}{\Psi \; u_{i}^{n}\frac{\partial}{\partial x_{i}}\left( {u_{j}^{n}\frac{\partial\varphi^{n}}{\partial x_{j}}} \right)\ {\Omega}}}}}\end{matrix}} & (19)\end{matrix}$

An embodiment of the present invention may include reformulatingequation (19) by integrating the last term of equation (19) by parts andemploying the divergence theorem to produce equation (20). In whichn_(i) is representative of the components of a normal vector {rightarrow over (n)}(Γ) that is perpendicular to a boundary Γ whichencapsulates integration space Ω.

$\begin{matrix}{{\int_{\Omega}^{\;}{\Psi \; \Psi^{T}\ {{\Omega \left( {\Phi^{n + 1} - \Phi^{n}} \right)}}}} = \begin{matrix}{{{- \Delta}\; t{\int_{\Omega}^{\;}{\sum\limits_{i}^{\;}{\Psi \; u_{i}^{n}\frac{\partial\varphi^{n}}{\partial x_{i}}\ {\Omega}}}}} -} \\{{\frac{\Delta \; t^{2}}{2}{\int_{\Omega}^{\;}{\sum\limits_{i,j}^{\;}{\frac{\partial\Psi}{\partial x_{i}}u_{i}^{n}u_{j}^{n}\frac{\partial\varphi^{n}}{\partial x_{j}}\ {\Omega}}}}} +} \\{\frac{\Delta \; t^{2}}{2}{\int_{\Gamma}^{\;}{\sum\limits_{i,j}^{\;}{\Psi \; n_{i}u_{i}^{n}u_{j}^{n}\frac{\partial\varphi^{n}}{\partial x_{j}}\ {\Gamma}}}}}\end{matrix}} & (20)\end{matrix}$

Equation (20) may be rearranged and rewritten in matrix form as inequation (21).

MΦ ^(n+1) =f=KΦ ^(n)   (21)

In which the matrices are defined in equations (22)-(23). Equation (23)describes is described in terms of a two-dimensional space. Anindividual skilled in the art would appreciate how equation (23) couldbe implemented as

$\begin{matrix}{M = {\int_{\Omega}^{\;}{\Psi \; \Psi^{T}\ {\Omega}}}} & (22) \\{K = {{\int_{\Omega}^{\;}{\Psi \; \Psi^{T}\ {\Omega}}} - {\Delta \; t{\int_{\Omega}^{\;}{\Psi \; B^{T}\ {\Omega}}}} - {\frac{\Delta \; t^{2}}{2}{\int_{\Omega}^{\;}{{BB}^{T}\ {\Omega}}}} + {\frac{\Delta \; t^{2}}{2}{\int_{\Gamma}^{\;}{{\Psi \left( {{n_{1}u_{1}^{n}} + {n_{2}u_{2}^{n}}} \right)}B^{T}\ {\Gamma}}}}}} & (23)\end{matrix}$

The following integrals in equations (24)-(26) can be used to describethe matrix B in terms of integral equations

$\begin{matrix}{{\int_{\Omega}^{\;}{{\Psi \left( {{u_{1}^{n}\frac{\partial\varphi^{n}}{\partial x_{1}}} + {u_{2}^{n}\frac{\partial\varphi^{n}}{\partial x_{2}}}} \right)}\ {\Omega}}} = {\int_{\Omega}^{\;}{\Psi \; B^{T}\ {\Omega}\; \Phi^{n}}}} & (24) \\{{\int_{\Gamma}^{\;}{{\Psi \left( {{n_{1}u_{1}^{n}} + {n_{2}u_{2}^{n}}} \right)}\left( \ {{u_{1}^{n}\frac{\partial\varphi^{n}}{\partial x_{1}}} + {u_{2}^{n}\frac{\partial\varphi^{n}}{\partial x_{2}}}} \right){\Gamma}}} = {\int_{\Gamma}^{\;}{{\Psi \left( {{n_{1}u_{1}^{n}} + {n_{2}u_{2}^{n}}} \right)}B^{T}\ {\Gamma}\; \Phi^{n}}}} & (25) \\{{{\int_{\Omega}^{\;}{\left( {{u_{1}^{n}\frac{\partial\Psi}{\partial x_{1}}} + {u_{2}^{n}\frac{\partial\Psi}{\partial x_{2}}}} \right)\left( {{u_{1}^{n}\frac{\partial\varphi^{n}}{\partial x_{1}}} + {u_{2}^{n}\frac{\partial\varphi^{n}}{\partial x_{2}}}} \right)\ {\Omega}}} = {\int_{\Omega}^{\;}{{BB}^{T}\ {\Omega}\; \Phi^{n}}}}{B = {\sum\limits_{i}^{\;}{u_{i}^{n}\frac{\partial\Psi}{\partial x_{i}}}}}} & (26)\end{matrix}$

Reinitializing the Level Set Function

An embodiment of the present invention uses the meshfree method tosimulate the motion of the level set function by using the level setevolution equation. The new level set function produced by the presentinvention m no longer be a signed distance function such that |∇φ|≠1.One method of addressing this issue is to re-initialize the level setfunction frequently such that, |∇φ|=1. Conventional routines forreinitializing a distance function enforce φ=0 at the interface andreset φ at all points close to the interface, which requires finding theinterface and significantly increases the resources required to simulatemany systems. One method for solving this problem is to integrate areinitialization equation such as the one shown in equation (27).

φ_(τ) +s(φ₀)|∇φ−1|=0   (27)

Equation (27) may be integrated for a short period of artificial time τ,where φ shares the same zero level set with the current level setfunction φ₀, τ is an artificial time period, s(φ₀) is the smoothed signfunction and φ₀ is the current level set to be reinitialized. Thecurrent level set function φ₀ may be far from the signed distancefunction.

The reinitialization of the level set function may be achieved bysolving the following partial differential equation (27). Until a steadystate form of the level set function φ is reached everywhere or atpoints near the interface. Equation (27) may be written in a steadystate form as in equation (28). The terms used in equation (28) may befound in equations (29) and equation (30).

$\begin{matrix}{{\frac{\partial\varphi}{\partial\tau} + {\overset{->}{c} \cdot {\nabla\; \varphi}}} = {{\frac{\partial\varphi}{\partial\tau} + {\sum\limits_{i}^{\;}{c_{i}\frac{\partial\varphi}{\partial x_{i}}}}} = {s\left( \varphi_{0} \right)}}} & (28) \\{{\overset{->}{c} = {{s\left( \varphi_{0} \right)}\frac{\nabla\varphi}{{\nabla\varphi}}}},{{s\left( \varphi_{0} \right)} = {s_{0} = \frac{\varphi_{0}}{\sqrt{\varphi_{0}^{2} + \left( {{{\nabla\varphi_{0}}}ɛ} \right)^{2}}}}}} & (29) \\{{\varphi_{0}\left( {\overset{->}{x},{\tau = 0}} \right)} = {\varphi\left( {\overset{->}{x},t} \right)}} & (30)\end{matrix}$

Equation (28) may be rewritten as equation (31).

$\begin{matrix}{\frac{\partial\varphi}{\partial\tau} = {{s\left( \varphi_{0} \right)} - {\sum\limits_{i}^{\;}{c_{i}\frac{\partial\varphi}{\partial x_{i}}}}}} & (31)\end{matrix}$

Equation (32) shows how the second temporal derivative of the level setfunction φ may be written using equation (31).

$\begin{matrix}\begin{matrix}{\frac{\partial^{2}\varphi}{\partial\tau^{2}} = {\frac{\partial\;}{\partial\tau}\left( {{s\left( \varphi_{0} \right)} - {\sum\limits_{i}^{\;}{c_{i}\frac{\partial\varphi}{\partial x_{i}}}}} \right)}} \\{= {- {\sum\limits_{i}^{\;}{c_{i}\frac{\partial}{\partial\tau}\left( \frac{\partial\varphi}{\partial x_{i}} \right)}}}} \\{= {- {\sum\limits_{i}^{\;}{c_{i}\frac{\partial}{\partial x_{i}}\left( \frac{\partial\varphi}{\partial\tau} \right)}}}} \\{= {- {\sum\limits_{i}^{\;}{c_{i}\frac{\partial}{\partial x_{i}}\left( {{s\left( \varphi_{0} \right)} - {\sum\limits_{j}^{\;}{c_{j}\frac{\partial\varphi}{\partial x_{j}}}}} \right)}}}} \\{= {\sum\limits_{i,j}^{\;}{c_{i}\frac{\partial}{\partial x_{i}}\left( {c_{j}\frac{\partial\varphi}{\partial x_{j}}} \right)}}}\end{matrix} & (32)\end{matrix}$

Equations (31) and (32) may be substituted into equation (13) to formequation (33). Equation (33) describes the temporal evolution of thelevel set function from a first point in time to a second point in timein terms of spatial derivative of the level set function at a firstpoint in time.

$\begin{matrix}\begin{matrix}{{\varphi^{n + 1} - \varphi^{n}} = {{\Delta \; \tau \frac{\partial\varphi^{n}}{\partial\tau}} + {\frac{\Delta \; \tau^{2}}{2}\frac{\partial\varphi^{n^{2}}}{\partial\tau^{2}}}}} \\{= {{\sum\limits_{i,j}^{\;}{\Delta \; {\tau \left( {s_{0} - {c_{i}^{n}\frac{\partial\varphi^{n}}{\partial x_{i}}}} \right)}}} + {\frac{\Delta \; \tau^{2}}{2}c_{i}^{n}\frac{\partial\;}{\partial x_{i}}\left( {c_{j}^{n}\frac{\partial\varphi^{n}}{\partial x_{j}}} \right)}}}\end{matrix} & (33)\end{matrix}$

The relations in equations (34) may be derived from equation (29).

$\begin{matrix}{{{\sum\limits_{i}^{\;}{c_{i}\frac{\partial\varphi}{\partial x_{i}}}} = {{s_{0}{\frac{\nabla\varphi}{{\nabla\varphi}} \cdot {\nabla\varphi}}} = {s_{0}{{\nabla\varphi}}}}}{and}{{\sum\limits_{i}^{\;}{c_{i}s_{0}{{\nabla\varphi}}}} = {{s_{0}\frac{\nabla\varphi}{{\nabla\varphi}}s_{0}{{\nabla\varphi}}} = {\nabla\varphi}}}} & (34)\end{matrix}$

Substituting the relations from Equation (34) into equation (33) givesus equation (35)

$\begin{matrix}\begin{matrix}{{\varphi^{n + 1} - \varphi^{n}} = {{\sum\limits_{i,j}^{\;}{\Delta \; {\tau \left( {s_{0} - {c_{i}^{n}\frac{\partial\varphi^{n}}{\partial x_{i}}}} \right)}}} + {\frac{\Delta \; \tau^{2}}{2}c_{i}^{n}\frac{\partial}{\partial x_{i}}\left( {c_{j}^{n}\frac{\partial\varphi^{n}}{\partial x_{j}}} \right)}}} \\{= {{\sum\limits_{i}^{\;}{\Delta \; {\tau \left( {s_{0} - {s_{0}{{\nabla\varphi^{n}}}}} \right)}}} + {\frac{\Delta \; \tau^{2}}{2}c_{i}^{n}\frac{\partial}{\partial x_{i}}\left( {s_{0}{{\nabla\varphi^{n}}}} \right)}}} \\{= {{\sum\limits_{i}^{\;}{\Delta \; {\tau \left( {s_{0} - {s_{0}{{\nabla\varphi^{n}}}}} \right)}}} + {\frac{\Delta \; \tau^{2}}{2}\begin{Bmatrix}{{\frac{\partial\;}{\partial x_{i}}\left( {c_{i}^{n}s_{0}{{\nabla\varphi^{n}}}} \right)} -} \\{s_{0}\frac{\partial c_{i}^{n}}{\partial x_{i}}{{\nabla\varphi^{n}}}}\end{Bmatrix}}}} \\{= {{\sum\limits_{i}^{\;}{\Delta \; \tau \left( {s_{0} - {s_{0}{{\nabla\varphi^{n}}}}} \right)}} + {\frac{\Delta \; \tau^{2}}{2}\left\{ {\frac{\partial{\nabla\varphi^{n}}}{\partial x_{i}} - {s_{0}\frac{\partial c_{i}^{n}}{\partial x_{i}}{{\nabla\varphi^{n}}}}} \right\}}}} \\{= {{\sum\limits_{i}^{\;}{\Delta \; \tau \left( {s_{0} - {s_{0}{{\nabla\varphi^{n}}}}} \right)}} - {\frac{\Delta \; \tau^{2}}{2}s_{0}\frac{\partial c_{i}^{n}}{\partial x_{i}}{{\nabla\varphi^{n}}}} +}} \\{{\frac{\Delta \; \tau^{2}}{2}\frac{\partial{\nabla\varphi^{n}}}{\partial x_{i}}}}\end{matrix} & (35)\end{matrix}$

As before the level set function may be rewritten as in equation (18) interms of a test function Ψ. Equation (18) may than be substituted intoequation (35) and integrated over the solution space in terms of thetest function Ψ, as shown in equation (36).

$\begin{matrix}{{\int_{\Omega}^{\;}{\Psi \; \Psi^{T}\ {{\Omega \left( {\Phi^{n + 1} - \Phi^{n}} \right)}}}} = {{\Delta \; \tau {\int_{\Omega}^{\;}{\Psi \; {s_{0}\left( {1 - {{\nabla\varphi^{n}}}} \right)}\ {\Omega}}}} - {\frac{\Delta \; \tau^{2}}{2}{\int_{\Omega}^{\;}{\Psi \; s_{0}\frac{\partial c_{i}^{n}}{\partial x_{i}}{{\nabla\varphi^{n}}}\ {\Omega}}}} - {\frac{\Delta \; \tau^{2}}{2}{\int_{\Omega}^{\;}{\frac{\partial\Psi}{\partial x_{i}}\frac{\partial\varphi^{n}}{\partial x_{i}}\ {\Omega}}}} + {\frac{\Delta \; \tau^{2}}{2}{\int_{\Gamma}^{\;}{\Psi \frac{\partial\varphi^{n}}{\partial x_{i}}n_{i}\ {\Gamma}}}}}} & (36)\end{matrix}$

As with equation (19) the surface Γ is the boundary that encapsulatesthe integration space Ω. The spatial derivate of the vector {right arrowover (c)}, δc_(i)i/δx_(i), is also a second-order spatial derivative ofthe level set function φ. If linear interpolation functions are used todescribe the level set function than δc_(i)/δx_(i) can be omitted.

Equation (36) is a nonlinear hyperbolic equation. The level set functionφ may be reinitialized so that |∇φ|=1, at least at points in thesimulation space Ω near the interface S. In an embodiment of the presentinvention, it may only be necessary for the level set function to behaveas an accurate distance function at points near the interface S. It mayonly be necessary to integrate equation (36) over τ=0 to t=aΔx, in whichΔx=mesh size and a=0.8. In an embodiment of the present invention, theboundary integral of equation (36) is not significant and may beignored.

Equation (37) is a simplified version of equation (36) in which a linearinterpolation of the function is assumed, and the boundary integral overthe boundary Γ is ignored.

$\begin{matrix}{{\int_{\Omega}^{\;}{\Psi \; \Psi^{T}\ {\Omega}\; \Phi^{m + 1}}} = {{\Delta \; \tau \; {\int_{\Omega}^{\;}{\Psi \; {s_{0}\left( {1 - {{\nabla\; \varphi^{m}}}} \right)}\ {\Omega}}}} + {\begin{bmatrix}{{\int_{\Omega}^{\;}{\Psi \; \Psi^{T}\ {\Omega}}} -} \\{\frac{{\Delta \; \tau^{2}}\;}{2}{\int_{\Omega}^{\;}{\frac{\partial\Psi}{\partial x_{i}}\frac{\partial\Psi^{T}}{\partial x_{i}}\ {\Omega}}}}\end{bmatrix}\Phi^{m}}}} & (37)\end{matrix}$

Inclusion of Boundary Terms

Equations (20) and (36) include boundary integrals along the boundary Γ.As noted above an embodiment of the present invention may includeignoring the effect of the boundary Γ. An alterative embodiment of thepresent invention may include integrating along the boundary Γ, whichwould improve the accuracy of solutions found using the presentinvention.

Equations (38) and (39) are the boundary integrals in equations (20) and(36) that were ignored in embodiments of the present invention describedabove.

$\begin{matrix}{\int_{\Gamma}^{\;}{\sum\limits_{i,j}^{\;}{\Psi \; n_{i}u_{i}^{n}u_{j}^{n}\frac{\partial\varphi^{n}}{\partial x_{j}}\ {\Gamma}}}} & (38) \\{\int_{\Gamma}^{\;}{\sum\limits_{i}^{\;}{\Psi \; \frac{\partial\varphi^{n}}{\partial x_{i}}n_{i}\ {\Gamma}}}} & (39)\end{matrix}$

If the problem domain Ω is two dimensional than the integrals (38) and(39) are line integrals. While, if the problem domain Ω is threedimensional than the integrals (38) and (39) are boundary integrals. Anindividual skilled in the art would appreciate how to adapt a twodimensional solution to a three dimensional solution, without goingbeyond the scope and spirit of the present invention as described in theclaims.

Since the nodal values of Ω^(n) at every time step is known to getφ^(n+1), Eqn. (3.7) and (3.14); therefore, φ^(n) in general can be foundin terms of RK shape function, the derivatives can be easily obtainedusing the derivatives of the RK shape functions where NP is the numberof discrete nodes.

$\begin{matrix}{\varphi^{h} = {\sum\limits_{I = 1}^{NP}{{\Psi_{I}\left( \overset{->}{x} \right)}\varphi_{I}}}} & (40) \\{{\frac{\partial\varphi^{h}}{\partial x_{j}} = {\sum\limits_{I = 1}^{NP}{\frac{\partial{\Psi_{I}\left( \overset{->}{x} \right)}}{\partial x_{j}}\varphi_{I}}}}{{j = x},y,z}} & (41) \\{{\frac{\partial\varphi^{h}}{\partial x_{j}} = {{\sum\limits_{I = 1}^{NP}{\frac{\partial{\Psi_{I}\left( \overset{->}{x} \right)}}{\partial x_{j}}\varphi_{I}}} + {\Psi_{I}\frac{\partial{\varphi_{I}\left( \overset{->}{x} \right)}}{\partial x_{j}}}}}{{j = x},y,z}} & (41)\end{matrix}$

FIG. 7 is an illustration of a boundary. The vector {right arrow over(n)} is normal to the boundary as shown in FIG. 7. The boundary may bedescribed using a series of nodes. The components of the vector {rightarrow over (n)} may be described in terms of its components as stated inequation (42). These components may be calculated from the position ofthe nodes as described in equations (43)-(45).

$\begin{matrix}{\overset{->}{n} = {n_{x,i} + n_{y,i}}} & (42) \\{L_{12} = \sqrt{\left( {x_{1} - x_{2}} \right)^{2} + \left( {y_{1} - y_{2}} \right)^{2}}} & (43) \\{n_{x,12} = {\left( {y_{2} - y_{1}} \right)/L_{12}}} & (44) \\{n_{y,12} = {\left( {x_{1} - x_{2}} \right)/L_{12}}} & (45)\end{matrix}$

An embodiment of the present invention is distinguished from the priorart at least by making use of an integral along the boundary whenreinitializing the boundary, as described in nonlinear equation (36).Prior art methods include the direct evaluation method that seeks outthe zero level set Γ and recalculates the signed-distance function φeverywhere in the solution domain Ω. However, these prior art methodsintroduce considerable complication to the evolution andreinitialization procedures and require significant computationalresources. The present invention lowers the resource requirements bymaking use of consistent integral equations. Computational resources arefurther reduced formulating the problem in terms of linear matrixequations. The present invention allows for these integral equations tobe solved effectively and consistently using Reproducing Kernel shapefunctions.

Calculating the interface integrals (38) and (39) includes calculatingthe spatial derivative of the level set function. Calculating spatialderivatives of the level set function includes calculating spatialderivatives of the shape functions. When the shape function is a linearfunction the spatial derivate is a constant. An embodiment of thepresent invention may benefit from shape functions in which the spatialderivate of the shape function is not constant over the integrationspace. The shape function may take the form of a polynomial or someother non-linear function. In a preferred embodiment of the presentinvention, a second order shape function is used. As the dimension ofthe problem space increases, the advantage of a shape function that isnot a first order approximation increases. The advantage of a higherorder shape function also increases as the complexity of the interfaceincrease.

Numerical Results

FIGS. 2A-2C show the evolution of a one-dimensional level set functionover a series of time steps. FIG. 2A is an illustration of aone-dimensional level set function that is described using two-nodelinear finite elements. The time step in the evolution equation is 0.072seconds. FIG. 2A show an initial level set function that at time t=0 isvery smooth but after several latter time steps is no longer smooth.This is an artifact of the simulation and is not representative of thesystem be simulated.

FIG. 2B is an illustration of the evolution of a one-dimensional levelset function like the one simulated in FIG. 2A. In FIG. 2B the level setfunction is described using three-node linear finite elements. The timestep in this case is 0.032 seconds. As with the previous case, unstablesolutions become clearly visible as time progresses.

FIG. 2C is an illustration of the evolution of a one-dimensional levelset function in which the present invention was used to simulate themotion of the level set function. As FIG. 2C illustrates the results ofthe present invention are noticeably better results than the resultsusing the finite element shown in FIGS. 2A-B.

FIGS. 3A-3E are illustrations of the results of an embodiment of thepresent invention used to simulate the evolution of a two-dimensionallevel set function. In this simulation u₁=u₂=1.0 and Δt=0.01 and: FIG.3A is at time t=0.0; FIG. 3B is at time t=0.02; FIG. 3C is at timet=0.04; FIG. 3D is at time t=0.06; and FIG. 3E is at time t=0.08;

FIG. 4A is an illustration of the evolution of a one dimensional levelset function φ at several different time steps as used in an embodimentof the present invention. FIG. 4B is an illustration of a spatialderivative of the level set function as shown in FIG. 4A. FIGS. 4A-B areillustrations of the evolution of one dimensional level set function inwhich the level set is not reinitialized and the boundary integral isignored.

FIG. 5A is an illustration of the evolution of a one dimensional levelset function φ at several different time steps in which the level setfunction is reinitialized as used in an embodiment of the presentinvention. FIG. 5B is an illustration of a spatial derivative of thelevel set function as shown in FIG. 5A. FIGS. 5A-B are illustrations ofthe evolution of one dimensional level set function in which the levelset function is reinitialized and the boundary integral is ignored.

FIGS. 6A-F are illustrations of the evolution of a two dimensional levelset function φ at several different time steps in which the level setfunction is reinitialized as used in an embodiment of the presentinvention. FIG. 6A is at time step t=0. FIG. 6B is at time step t=0.25.FIG. 6C is at time step t=0.50. FIG. 6D is at time step t=0.75. FIG. 6Eis at time step t=1.00. FIG. 6F is at time step t=1.25.

FIG. 8A is an illustration of the evolution of a one dimensional levelset function φ at several different time steps in which the level setfunction is reinitialized and the boundary integral is not ignored asused in a preferred embodiment of the present invention. FIG. 8B is anillustration of a spatial derivative of the level set function as shownin FIG. 8A.

Applications

An embodiment of the present invention may be used to describe aninterface. The interface may be between a first material and a secondmaterial. The interface may be between a first fluid and a second fluid.The first fluid may be ink. The second fluid may be air. An embodimentof the present invention may also be used to describe the motion of theinterface. The simulation of the motion of the interface may be used ina simulation of an ink jet print head. An embodiment of the presentinvention may be used in the process of designing an ink jet head. Anembodiment of the present invention may be used in the process ofdesigning an object in which the motion of the interface is a matter ofinterest. An embodiment of the present invention may be an object thatwas designed using the method described above as a guide. An individualskilled in the art would appreciate that the present invention may beused to design other products.

System

Having described the details of the invention, an exemplary system 1000,which may be used to implement one or more aspects of the presentinvention will now be described with reference to FIG. 9. As illustratedin FIG. 9, the system includes a central processing unit (CPU) 1001 thatprovides computing resources and controls the computer. The CPU 1001 maybe implemented with a microprocessor or the like, and may also include agraphics processor and/or a floating point coprocessor for mathematicalcomputations. The system 1000 may also include system memory 1002, whichmay be in the form of random-access memory (RAM) and read-only memory(ROM).

A number of controllers and peripheral devices may also be provided, asshown in FIG. 9. An input controller 1003 represents an interface tovarious input device(s) 1004, such as a keyboard, mouse, or stylus.There may also be a scanner controller 1005, which communicates with ascanner 1006. The system 1000 may also include a storage controller 1007for interfacing with one or more storage devices 1008 each of whichincludes a storage medium such as magnetic tape or disk, or an opticalmedium that might be used to record programs of instructions foroperating systems, utilities and applications which may includeembodiments of programs that implement various aspects of the presentinvention. Storage device(s) 1008 may also be used to store processeddata or data to be processed in accordance with the invention. Thesystem 1000 may also include a display controller 1009 for providing aninterface to a display device 1011, which may be a cathode ray tube(CRT), or a thin film transistor (TFT) display. The system 1000 may alsoinclude a printer controller 1012 for communicating with a printer 1013.A communications controller 1014 may interface with one or morecommunication devices 1015 which enables the system 1000 to connect toremote devices through any of a variety of networks including theInternet, a local area network (LAN), a wide area network (WAN), orthrough any suitable electromagnetic carrier signals including infraredsignals.

In the illustrated system, all major system components may connect to abus 1016, which may represent more than one physical bus. However,various system components may or may not be in physical proximity to oneanother. For example, input data and/or output data may be remotelytransmitted from one physical location to another. In addition, programsthat implement various aspects of this invention may be accessed from aremote location (e.g., a server) over a network. Such data and/orprograms may be conveyed through any of a variety of machine-readablemedium including magnetic tape or disk or optical disc, or atransmitter, receiver pair.

The present invention may be conveniently implemented with software.However, alternative implementations are certainly possible, including ahardware implementation or a software/hardware implementation. Anyhardware-implemented functions may be realized using ASIC(s), digitalsignal processing circuitry, or the like. Accordingly, the “means” termsin the claims are intended to cover both software and hardwareimplementations. Similarly, the term “computer-readable medium” as usedherein includes software and or hardware having a program ofinstructions embodied thereon, or a combination thereof. With theseimplementation alternatives in mind, it is to be understood that thefigures and accompanying description provide the functional informationone skilled in the art would require to write program code (i.e.,software) or to fabricate circuits (i.e., hardware) to perform theprocessing required.

In accordance with further aspects of the invention, any of theabove-described methods or steps thereof may be embodied in a program ofinstructions (e.g., software), which may be stored on, or conveyed to, acomputer or other processor-controlled device for execution on acomputer readable medium. Alternatively, any of the methods or stepsthereof may be implemented using functionally equivalent hardware (e.g.,application specific integrated circuit (ASIC), digital signalprocessing circuitry, etc.) or a combination of software and hardware.

While the invention has been described in conjunction with severalspecific embodiments, it is evident to those skilled in the art thatmany further alternatives, modifications, and variations will beapparent in light of the foregoing description. Thus, the inventiondescribed herein is intended to embrace all such alternatives,modifications, applications and variations as may fall within the spiritand scope of the appended claims.

1. A method for simulating the motion of an interface moving through asimulation space comprising: simulating the interface using a level setfunction to describe a position and shape of the interface in thesimulation space at a first point in time describing the level setfunction at the first point time using a meshfree method; describing amotion of the interface from the first point in time to a second pointtime using a level set evolution method; and finding an approximatesolution to the level set evolution method using the meshfree method todescribe the level set function at the second point in time.
 2. Themethod of claim 1, wherein the meshfree method is a Reproducing KernelParticle Method.
 3. The method of claim 1, wherein the interface isbetween a first fluid and a second fluid.
 4. The method of claim 3,wherein the first fluid is ink, and the second fluid is air.
 5. Acomputer-readable medium encoded with instructions to a computer forperforming the method of claim
 1. 6. A system including a memory moduleand instructions for performing the method of claim
 1. 7. The method ofclaim 1, wherein the meshfree method approximates the level set functionusing a first family of functions that define the smoothness of theapproximation of the level set function, and each function in the firstfamily of functions has a compact support, and the size of the compactsupport is a function of a maximum of a set of k-th ordered statisticsof a plurality of distances between sets of distinct nodes in thesimulation space.
 8. The method of claim 7, wherein k is selected fromthe group consisting of 2, 3, 4, and
 5. 9. The method of claim 7,wherein the size of the compact support is also a function of an orderof an enrichment function.
 10. The method of claim 1, wherein the levelset evolution method is in described in terms of spatial derivates ofthe level set function; the spatial derivatives of the level setfunction are discretized, and a mesh free method is used to describeeach of the discretized terms of the spatial derivatives of the levelset function.
 11. The method of claim 10, wherein the level setevolution method includes an integral along a boundary surrounding thesimulation space, including a spatial derivate of the level setfunction, a test function, and components of a vector normal to theboundary.
 12. The method of claim 1, further comprising the designing aphysical object based on the results of the method of simulating themotion of an interface.
 13. The method of claim 12, wherein the physicalobject is an ink jet head.
 14. The method of claim 7 wherein eachfunction in the family of function is positive in the area bounded bythe size of the compact support.
 15. The method of claim 1, wherein themeshfree method approximates the level set function using a first familyof functions that define the smoothness of the approximation of thelevel set function, and each function in the first family of functionshas a compact support, and the size of the compact support is a functionof a maximum of a set of k-th ordered statistics of a plurality of powersums of differences between components that define the positions of setsof distinct nodes in the simulation space.
 16. The physical objectdesigned using the method of claim
 12. 17. The physical object of claim16, wherein the physical object is an ink jet head.